package org.chorus.utils;

// *********************************************************************************************
/** Simple Kernel Density Estimator that uses one Gaussian kernel per observed data value.
 *  This class is adapted from the Weka software package (weka.sourceforge.net).
 *  @author Len Trig  (trigg@cs.waikato.ac.nz) with modifications by Jo Wood, giCentre.
 *  @version 1.0, 24th May, 2011.
 */ 
// *********************************************************************************************

/* This file is part of the Chorus river bed analysis software package. Chorus is free software:
 * you can redistribute it and/or modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 * 
 * Chorus is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
 * See the GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along with this
 * source code (see COPYING included with this source code). If not, see 
 * http://www.gnu.org/licenses/.
 */

public class KernelEstimator 
{
	// -------------------------------- Object variables --------------------------------

	private static double SMALL = 1e-6; 	// The small deviation allowed in double comparisons.
	private static double MAX_ERROR = 0.01;	// Maximum percentage error permitted in probability calculations.
	private static final double SQRTH  =  7.07106781186547524401E-1;
	private static final double MAXLOG =  7.09782712893383996732E2;
	
	private double[] values;		// Original data values from which to base estimates.
	private int numValues;			// Current number of values in estimator.
	private double precision;		// The precision of data values.
	private double bandwidth;		// Kernel bandwidth (smoothing).

	// ---------------------------------- Constructors ----------------------------------

	/** Creates an estimator with the given bandwidth. The bandwidth corresponds to the 
	 *  standard deviation of the Gaussian kernel used to estimate probabilities. Larger
	 *  bandwidths produce smoother estimates. This version of the constructor assumes a
	 *  precision of 1 for all data values (i.e. all are rounded to their nearest integer
	 *  category).
	 *  @param bandwidth Bandwidth of the estimator in the same units as data values.
	 */
	public KernelEstimator(double bandwidth)
	{
		this(bandwidth,1);
	}

	/** Creates an estimator with the given bandwidth and precision. The bandwidth 
	 *  corresponds to the standard deviation of the Gaussian kernel used to estimate 
	 *  probabilities. Larger bandwidths produce smoother estimates. The precision value
	 *  controls the degree of rounding applied to the data. For example, a precision of
	 *  0.1, treats values in the interval (0.25,0.35] as 0.3.
	 *  @param bandwidth Bandwidth of the estimator in the same units as data values.
	 *  @param precision The precision of the data values to be smoothed.
	 */
	public KernelEstimator(double bandwidth, double precision)
	{
		this.bandwidth = bandwidth;
		values = new double [50];
		numValues = 0;
		this.precision = Math.max(precision, SMALL);
	}

	// ------------------------------------ Methods -------------------------------------


	/** Adds a new data value to the current estimator.
	 *  @param value New data value to add. 
	 */
	public void addValue(double value) 
	{
		double roundedValue = round(value);
		int insertIndex = findNearestValue(roundedValue);
		if ((numValues <= insertIndex) || (values[insertIndex] != roundedValue))
		{
			if (numValues < values.length) 
			{
				int left = numValues - insertIndex; 
				System.arraycopy(values, insertIndex, values, insertIndex + 1, left);
				values[insertIndex] = roundedValue;
				numValues++;
			} 
			else
			{
				double [] newValues = new double [values.length * 2];
				int left = numValues - insertIndex; 
				System.arraycopy(values, 0, newValues, 0, insertIndex);
				newValues[insertIndex] = roundedValue;
				System.arraycopy(values, insertIndex, newValues, insertIndex + 1, left);
				numValues++;
				values = newValues;
			}
		}
	}

	/** Reports the probability estimate for a given value. Given that probabilities of any 
	 *  discrete value is 0, this probability assumes binning to the precision identified in
	 *  the constructor.
	 *  @param value The value whose probability is to be estimated. 
	 *  @return Estimated probability of the supplied value
	 */
	public double getProbability(double value)
	{
		double delta = 0, sum = 0, currentProb = 0;
		double zLower = 0, zUpper = 0;
		if (numValues == 0) 
		{
			zLower = (value - (precision / 2)) / bandwidth;
			zUpper = (value + (precision / 2)) / bandwidth;
			return (normalProbability(zUpper) - normalProbability(zLower));
		}

		double weightSum = 0;
		int start = findNearestValue(value);
		for (int i=start; i<numValues; i++) 
		{
			delta = values[i] - value;
			zLower = (delta - (precision / 2)) / bandwidth;
			zUpper = (delta + (precision / 2)) / bandwidth;
			currentProb = (normalProbability(zUpper) - normalProbability(zLower));
			sum += currentProb;
			weightSum++;
			if (currentProb*(numValues-weightSum) < sum*MAX_ERROR) 
			{
				break;
			}
		}

		for (int i=start-1; i>=0; i--) 
		{
			delta = values[i] - value;
			zLower = (delta - (precision / 2)) / bandwidth;
			zUpper = (delta + (precision / 2)) / bandwidth;
			currentProb = (normalProbability(zUpper) - normalProbability(zLower));
			sum += currentProb;
			weightSum ++;
			if (currentProb*(numValues-weightSum) < sum*MAX_ERROR) 
			{
				break;
			}
		}
		return sum / numValues;
	}


	/** Provides a textual summary of this estimator 
	 *  @return Text representing this estimator.
	 */
	public String toString() 
	{
		String result = numValues+" Normal Kernels. \nbandwidth="+doubleToString(bandwidth,6,4)+" Precision="+precision;
		if (numValues == 0) 
		{
			result += "  \nMean = 0";
		} 
		else 
		{
			result += "  \nMeans =";
			for (int i=0; i<numValues; i++)
			{
				result += " " + values[i];
			}
		}
		return result + "\n";
	}

	/** Reports the precision of this kernel estimator. This is the width of the discrete categories around 
	 *  each value that informed the estimator. 
	 *  @return Precision of the estimator as defined in the constructor.
	 */
	public double getPrecision() 
	{
		return precision;
	}

	/** Reports the bandwidth of this estimator. This is equivalent to the standard deviation of
	 *  the Gaussian distribution used to smooth data values. The larger the value, the more smoothing applied.
	 * @return Bandwidth of the estimator.
	 */
	public double getBandwidth() 
	{
		return bandwidth;
	}


	// --------------------------------- Private methods --------------------------------

	/** Reports the area under the Normal (Gaussian) probability density function, integrated from minus 
	 *  infinity to <tt>x</tt> (assumes mean is zero, variance is one).
	 * <pre>
	 *                            x
	 *                             -
	 *                   1        | |          2
	 *  normal(x)  = ---------    |    exp( - t /2 ) dt
	 *               sqrt(2pi)  | |
	 *                           -
	 *                          -inf.
	 *
	 *             =  ( 1 + erf(z) ) / 2
	 *             =  erfc(z) / 2
	 * </pre>
	 * where <tt>z = x/sqrt(2)</tt>.
	 * Computation is via the functions <tt>errorFunction</tt> and <tt>errorFunctionComplement</tt>.
	 * @param a the z-value
	 * @return the probability of the z value according to the normal pdf
	 */
	static double normalProbability(double a) 
	{ 
		double x, y, z;

		x = a * SQRTH;
		z = Math.abs(x);

		if (z<SQRTH)
		{
			y = 0.5 + 0.5*errorFunction(x);
		}
		else
		{
			y = 0.5*errorFunctionComplemented(z);
			if (x>0)
			{
				y = 1.0 - y;
			}
		} 
		return y;
	}

	/** Reports the error function of the normal distribution. The integral is
	 * <pre>
	 *                           x 
	 *                            -
	 *                 2         | |          2
	 *   erf(x)  =  --------     |    exp( - t  ) dt.
	 *              sqrt(pi)   | |
	 *                          -
	 *                           0
	 * </pre>
	 * <b>Implementation:</b>
	 * For <tt>0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2)</tt>; otherwise
	 * <tt>erf(x) = 1 - erfc(x)</tt>.
	 * <p>
	 * Code adapted from the <a href="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
	 * Java 2D Graph Package 2.4</a>, which in turn is a port from the
	 * <a href="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</a>
	 * Math Library (C).
	 * @param x the argument to the function.
	 */
	private static double errorFunction(double x)
	{ 
		double y, z;
		final double T[] =
		{
				9.60497373987051638749E0,
				9.00260197203842689217E1,
				2.23200534594684319226E3,
				7.00332514112805075473E3,
				5.55923013010394962768E4
		};
		final double U[] = 
		{
				3.35617141647503099647E1,
				5.21357949780152679795E2,
				4.59432382970980127987E3,
				2.26290000613890934246E4,
				4.92673942608635921086E4
		};

		if (Math.abs(x) > 1.0)
		{
			return(1.0 - errorFunctionComplemented(x));
		}

		z = x * x;
		y = x * polevl(z,T,4) / p1evl(z,U,5);
		return y;
	}

	/** Reports the complementary Error function of the normal distribution.
	 * <pre>
	 *  1 - erf(x) =
	 *
	 *                           inf. 
	 *                             -
	 *                  2         | |          2
	 *   erfc(x)  =  --------     |    exp( - t  ) dt
	 *               sqrt(pi)   | |
	 *                           -
	 *                            x
	 * </pre>
	 * <b>Implementation:</b>
	 * For small x, <tt>erfc(x) = 1 - erf(x)</tt>; otherwise rational approximations are computed.
	 * <p>
	 * Code adapted from the <a href="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
	 * Java 2D Graph Package 2.4</a>, which in turn is a port from the
	 * <a href="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</a>
	 * Math Library (C).
	 * @param a the argument to the function.
	 */
	private static double errorFunctionComplemented(double a)
	{ 
		double x,y,z,p,q;

		double P[] = 
		{
				2.46196981473530512524E-10,
				5.64189564831068821977E-1,
				7.46321056442269912687E0,
				4.86371970985681366614E1,
				1.96520832956077098242E2,
				5.26445194995477358631E2,
				9.34528527171957607540E2,
				1.02755188689515710272E3,
				5.57535335369399327526E2
		};
		double Q[] = 
		{
				1.32281951154744992508E1,
				8.67072140885989742329E1,
				3.54937778887819891062E2,
				9.75708501743205489753E2,
				1.82390916687909736289E3,
				2.24633760818710981792E3,
				1.65666309194161350182E3,
				5.57535340817727675546E2
		};

		double R[] = 
		{
				5.64189583547755073984E-1,
				1.27536670759978104416E0,
				5.01905042251180477414E0,
				6.16021097993053585195E0,
				7.40974269950448939160E0,
				2.97886665372100240670E0
		};
		double S[] = 
		{
				2.26052863220117276590E0,
				9.39603524938001434673E0,
				1.20489539808096656605E1,
				1.70814450747565897222E1,
				9.60896809063285878198E0,
				3.36907645100081516050E0
		};

		if (a<0.0)
		{
			x = -a;
		}
		else
		{
			x = a;
		}

		if(x<1.0)
		{
			return 1.0 - errorFunction(a);
		}

		z = -a * a;

		if(z<-MAXLOG) 
		{
			if (a<0)
			{
				return( 2.0 );
			}
			return( 0.0 );
		}

		z = Math.exp(z);

		if (x< 8.0)
		{
			p = polevl( x, P, 8 );
			q = p1evl( x, Q, 8 );
		}
		else
		{
			p = polevl( x, R, 5 );
			q = p1evl( x, S, 6 );
		}

		y = (z*p)/q;

		if (a<0)
		{
			y = 2.0-y;
		}

		if (y==0.0)
		{
			if (a<0)
			{
				return 2.0;
			}
			return(0.0);
		}
		return y;
	}

	/** Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
	 *  Evaluates polynomial when coefficient of N is 1.0. Otherwise same as <tt>polevl()</tt>.
	 *  <pre>
	 *                     2          N
	 * y  =  C  + C x + C x  +...+ C x
	 *        0    1     2          N
	 *
	 * Coefficients are stored in reverse order:
	 *
	 * coef[0] = C  , ..., coef[N] = C  .
	 *            N                   0
	 * </pre>
	 *  The function <tt>p1evl()</tt> assumes that <tt>coef[N] = 1.0</tt> and is
	 *  omitted from the array.  Its calling arguments are
	 *  otherwise the same as <tt>polevl()</tt>.
	 *  In the interest of speed, there are no checks for out of bounds arithmetic.
	 *  @param x argument to the polynomial.
	 *  @param coef the coefficients of the polynomial.
	 *  @param N the degree of the polynomial.
	 */
	private static double p1evl(double x, double coef[], int N)
	{
		double ans;
		ans = x + coef[0];

		for(int i=1; i<N; i++)
		{
			ans = ans*x+coef[i];
		}
		return ans;
	}

	/** Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
	 * <pre>
	 *                     2          N
	 * y  =  C  + C x + C x  +...+ C x
	 *        0    1     2          N
	 *
	 * Coefficients are stored in reverse order:
	 *
	 * coef[0] = C  , ..., coef[N] = C  .
	 *            N                   0
	 * </pre>
	 *  In the interest of speed, there are no checks for out of bounds arithmetic.
	 *  @param x argument to the polynomial.
	 *  @param coef the coefficients of the polynomial.
	 *  @param N the degree of the polynomial.
	 */
	private static double polevl(double x, double coef[], int N) 
	{
		double ans;
		ans = coef[0];

		for(int i=1; i<=N; i++)
		{
			ans = ans*x+coef[i];
		}
		return ans;
	}

	/** Rounds a double and converts it into a formatted decimal-justified String.
	 *  Trailing 0's are replaced with spaces.
	 *  @param value the double value
	 *  @param width the width of the string
	 *  @param afterDecimalPoint the number of digits after the decimal point
	 *  @return the double as a formatted string
	 */
	private static String doubleToString(double value, int width, int afterDecimalPoint)
	{
		String tempString = doubleToString(value, afterDecimalPoint);
		char[] result;
		int dotPosition;

		if ((afterDecimalPoint >= width) || (tempString.indexOf('E') != -1))
		{ 
			// Protects scientific notation
			return tempString;
		}

		// Initialize result
		result = new char[width];
		for (int i = 0; i < result.length; i++)
		{
			result[i] = ' ';
		}

		if (afterDecimalPoint > 0)
		{
			// Get position of decimal point and insert decimal point
			dotPosition = tempString.indexOf('.');
			if (dotPosition == -1) 
			{
				dotPosition = tempString.length();
			} 
			else
			{
				result[width - afterDecimalPoint - 1] = '.';
			}
		} 
		else
		{
			dotPosition = tempString.length();
		}

		int offset = width - afterDecimalPoint - dotPosition;
		if (afterDecimalPoint > 0)
		{
			offset--;
		}

		// Not enough room to decimal align within the supplied width
		if (offset < 0)
		{
			return tempString;
		}

		// Copy characters before decimal point
		for (int i = 0; i < dotPosition; i++)
		{
			result[offset + i] = tempString.charAt(i);
		}

		// Copy characters after decimal point
		for (int i = dotPosition + 1; i < tempString.length(); i++)
		{
			result[offset + i] = tempString.charAt(i);
		}
		return new String(result);
	}

	/** Rounds a double and converts it into String.
	 *  @param value the double value
	 *  @param afterDecimalPoint the (maximum) number of digits permitted
	 *  after the decimal point
	 *  @return the double as a formatted string
	 */
	private static String doubleToString(double value, int afterDecimalPoint)
	{
		StringBuffer stringBuffer;
		double temp;
		int dotPosition;
		long precisionValue;

		temp = value * Math.pow(10.0, afterDecimalPoint);
		if (Math.abs(temp) < Long.MAX_VALUE) 
		{
			precisionValue = (temp > 0) ? (long)(temp + 0.5) : -(long)(Math.abs(temp) + 0.5);
			if (precisionValue == 0)
			{
				stringBuffer = new StringBuffer(String.valueOf(0));
			} 
			else
			{
				stringBuffer = new StringBuffer(String.valueOf(precisionValue));
			}
			if (afterDecimalPoint == 0) 
			{
				return stringBuffer.toString();
			}
			dotPosition = stringBuffer.length() - afterDecimalPoint;
			while (((precisionValue < 0) && (dotPosition < 1)) || (dotPosition < 0)) 
			{
				if (precisionValue < 0) 
				{
					stringBuffer.insert(1, '0');
				} 
				else
				{
					stringBuffer.insert(0, '0');
				}
				dotPosition++;
			}
			stringBuffer.insert(dotPosition, '.');
			if ((precisionValue < 0) && (stringBuffer.charAt(1) == '.')) 
			{
				stringBuffer.insert(1, '0');
			}
			else if (stringBuffer.charAt(0) == '.') 
			{
				stringBuffer.insert(0, '0');
			}
			int currentPos = stringBuffer.length() - 1;
			while ((currentPos > dotPosition) && (stringBuffer.charAt(currentPos) == '0'))
			{
				stringBuffer.setCharAt(currentPos--, ' ');
			}
			if (stringBuffer.charAt(currentPos) == '.') 
			{
				stringBuffer.setCharAt(currentPos, ' ');
			}
			return stringBuffer.toString().trim();
		}
		return new String("" + value);
	}

	/** Locates the nearest data value using a binary search.
	 *  @param key Data value to locate
	 *  @return Index of the nearest data value
	 */
	private int findNearestValue(double key) 
	{
		int low = 0; 
		int high = numValues;
		int middle = 0;
		while (low < high)
		{
			middle = (low + high) / 2;
			double current = values[middle];
			if (current == key) 
			{
				return middle;
			}
			if (current > key) 
			{
				high = middle;
			}
			else if (current < key) 
			{
				low = middle + 1;
			}
		}
		return low;
	}

	/** Rounds a data value using the defined precision for this estimator
	 *  @param data the value to round
	 *  @return the rounded data value
	 */
	private double round(double data)
	{
		return Math.rint(data / precision) * precision;
	}
}
